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Abstract 

We study statistical properties of the irregular bursting arising in a class of neuronal models 
close to the transition from spiking to bursting. Prior to the transition to bursting, the systems 
in this class develop chaotic attractors, which generate irregular spiking. The chaotic spiking 
gives rise to irregular bursting. The duration of bursts near the transition can be very long. We 
describe the statistics of the number of spikes and the interspike interval distributions within 
one burst as functions of the distance from criticality. 

Bursting oscillations are ubiquitous in the experimental and modeling studies of excitable cell 
membranes. Many models generating bursting have been subject to intensive research due to 
their physiological significance and dynamical complexity (see ^ E3 El El and references 
therein). Under the variation of parameters even minimal 3D models of bursting neurons exhibit 
a rich variety of periodic and aperiodic dynamical patterns corresponding to different spiking and 
bursting regimes. The transitions between these patterns may contain complex dynamical struc- 
tures such as period-doubling cascades and deterministic chaos. In particular, it was shown that the 
transition from tonic spiking to bursting in a class of bursting neuron models, so-called square- wave 
bursters, contains windows of chaotic dynamics [H f2| HTH 115]. In view of the complex bifurcation 
structure of this class of problems, it is important to identify the universal features pertinent to 
different dynamical patterns and transitions between them. In the present Letter, we describe sta- 
tistical features of the irregular bursting arising in a class of neuronal models close to the transition 
from spiking to bursting. Prior to transition to bursting, the systems in this class develop chaotic 
attractors, which generate irregular spiking. The chaotic spiking gives rise to irregular bursting. 
The duration of bursts near the transition can be extremely long (see Figure We analyze the 
statistics of the number of spikes and the interspike intervals within one burst as functions of the 
distance from criticality. 

To describe our results, we use a three variable model of a bursting neuron introduced in 
[5]. The model is based on three nonlinear conductances: persistent sodium, InuP, the delayed 
rectifier, Ik, a slow potassium M-current, Irm, and a passive current, 1^. In spite of a number of 
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Figure 1: Periodic and aperiodic firing patterns generated by in the regimes close to the 

transition from spiking to bursting. The values of parameters are C = 1 (/Lf-Fcm~ 2 ) ; qna = 20, 
g K = 10, g L = 8 (mScm~ 2 ); E Na = 60, E K = -90, E L = -80 (mV); a m = -20, a n = -25, 
a w = —20 (mV); b m = 15, b n = 5, b w = 5; r n = 0.152, t w = 20 (rns -1 ), and / = hp A. 

simplifications, the model captures essential features of a class of more detailed physiological models 
and is representative for a class of square- wave bursters [5]. The latter are among most common 
models of bursting. The following system of three differential equations describes the dynamics of 
the membrane potential, v, and two gating variables n and w: 

Cv = F(v,n,w), (1) 

h = , (2) 

w = , (3) 

where F(v,n,w) = -5iVaP"i 00 (w)(w - E NaP ) - g K n{v - E K ) - -yw(v - E K ) - g L (v - E L ) + J; 
g s and E s , (s£ {NaP, K, L}) are the maximal conductance and the reversal potential of I s , s G 
{NaP, K, L} , respectively; and / is the applied current. The maximal conductance of Irm, 7, is 
viewed as a control parameter. The time constants r n and t w determine the rates of activation in 
the populations of K and KM channels. The steady state functions are defined by 

Soo{v) = — — — , se{m,n,w}. 

1 + exp 



The parameter values are given in the caption to Figure ^ 

The analysis of the bursting neuron models like ©-© uses a fast-slow decomposition [51 151 ITU]. 
Specifically, we note that the time constant t w presents the slowest time scale in the dynamics of 
Therefore, we view a = r" 1 > as a small parameter. In the limit as a — > 0, system 
0-© i s reduced to a 2D fast subsystem Q, (J2J), where w is viewed as a parameter. For values 
of w £ {waHiWhc) the fast subsystem has a family of stable periodic orbits, which is born at the 
Andronov-Hopf bifurcation at w = wah and terminates at the homoclinic bifurcation at w = whc- 
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Figure 2: The bifurcation diagram of the fast subsystem (JIJ and (J2J with a superimposed periodic 
trajectory induced by the slow equation (j3J). Curves S and C indicate the families of periodic orbits 
and the fixed points respectively. 

A trajectory of the full system ([! ) )- ([5 ]) approaches the surface foliated by the periodic orbits of the 
fast subsystem, S (Figure HJ). The evolution along S is governed by the slow equation: 

w = a (woo(v) - w) , (4) 

where the leading order approximation of v(t) is obtained from the fast equations and (p]). 
When the trajectory hits the boundary of S (i.e., when w > whc), it jumps down to the curve of 
stable fixed points C, the only attracting set of the fast subsystem for w > who- Then it evolves 
along C as shown in Figure [2] until it reaches the boundary of C at w = wsn, corresponding to a 
saddle-node (SN) bifurcation in the fast subsystem. This is followed by the jump back to S and 
the oscillations in the fast subsystem resume. This description accounts for one cycle of bursting 
oscillation. 

The description of bursting dynamics in a wide class of neuronal models, including (|T )l -l(3* ]> . can 
be reduced to a single difference equation for the slow variable. Indeed, since the state of the fast 
system is determined by the value of the slow variable w, it is sufficient to know how w changes 
after each cycle of oscillations of the fast subsystem and after the period of quiescence: 

w n+ i = P 7 (w n ) , n = 1, 2, (5) 

This idea underlies the method of reduction of a class of models of bursting neurons to one- 
dimensional maps proposed in [Tj- Similar map-based approaches were used in JTTJ for studying 
bursting patterns in the Belousov-Zhabotinskii reaction and in [Hj for analyzing complex oscillatory 
regimes in a compartmental model of the dopamine neuron. In |7|, we provided the analytical 
description for P 7 , which was used to account for various spiking and bursting patterns and transi- 
tions between in the class of models of excitable cells. For the model at hand, the one-dimensional 
map is shown in Figure 01 The bifurcation structure of the fast subsystem endows the map with 
distinct structure: it is a piecewise continuous map with the boundary layer, Jo, corresponding to 
the homoclinic bifurcation in the fast subsystem (Figure EJ). There are two intervals of continuity 
in the domain of definition of P 7 , I\ = I~ [J I and I2 = I + ■ The iterations of P 7 over I\ cor- 
respond to the changes of w after each spike of voltage and the definition of P 7 over I2 captures 
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Figure 3: The first return map for the slow variable w for the value of control parameter 7 = 2.9375 
close to the critical value 7°. 

the mechanism of return to spiking after the period of quiescence (Figure I3J) • In 1% , Py « wsn 
is almost constant and can be well approximated by a linear function with a small negative slope 
(see Remark 5.3(b) in [7]). By construction, P 7 has a unique fixed point w 7 G I\. For small values 
of 7, Wy is globally attracting. The trajectory of (JSJ) stays in a small neighborhood of u) 7 G /1 
after several iterations of P 7 . Therefore, the value of w in JIJ and (j2J changes little and the fast 
subsystem generates periodic spiking. The bifurcation scenario described in [7] for a model of a 
bursting neuron applies to a wide class problems including ([T )) -(|3* ]) . According to this scenario, for 
increasing values of 7 the fixed point w) 7 loses stability through a period-doubling (PD) bifurcation, 
giving rise to a stable period 2 orbit. In the continuous system Q-©, this corresponds to the 
appearance of the stable periodic solution formed by clustered pairs of spikes, so-called doublets. 
Under further increase of 7, the numerical experiments show other PD bifurcations of periodic 
orbits and, eventually, the dynamics of the discrete system © becomes chaotic. This scenario fits 
well with the numerical results reported for related models pfl El [JJ . While the maximum of P 7 
over Ii, P 7 = max we i 1 P 1 (w) remains less than whc, Pyih) C 1%, and the trajectories of © are 
trapped in I\. This means that, in this regime, the continuous system exhibits (possibly chaotic) 
spiking. The transition to bursting takes place at 7 = r y c : 



For 7 > 7 C , trajectories of J5J) may leave I\. However, for values of 7 just above 7 C , the window 
of escape, J 7 , is very small (see Figure 01 • Therefore, a trajectory of © with large probability 
spends a long time in I\ before escaping to Ii- If, in addition, the dynamics of © for 7 = 7 C has 
mixing property, the transition to bursting lies through regimes of chaotic bursting with very long 
intervals of spiking appearing with large probability. Below we study the statistics of the number 
of spikes and the interspike intervals within one burst for near the transition to bursting. 

Let I = [wq,whc] j wo = ]im w -> WHC -o Py(w) (see Figure 0} and assume that at the critical 
value of the control parameter 7 = 7 C , map P 7 c : I 1— ► I has an invariant probability measure 
H absolutely continuous to the Lebesgue measure on /. In addition, we assume that P 7 c has the 
mixing property: 
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Figure 4: The mean values of the number of spikes, iV 7 , are fitted with y = a (7 — 7 C ) 2 , (7 < 7) 
and y = 6(7 - 7 C ) _1 , (7 > 7), where a « 25.42, b 0.52, 7 C « 2.93419, and 7 2.934343. 

for all measurable sets A,B <Z I ^3]. Let 7 — 7 C > be sufficiently small. Denote 5 = /i(J 7 ), where 
J 7 = {u; G 7 : P~/(w) > whc} (Figure|3|). Then the expected value (with respect to the probability 
measure /i) of the number of spikes within one burst is given by 

00 

iV 7 = ^(J 7 ) + ^£;/i fc , (8) 

k=2 



where fi k = pi (P~ k+l ( J 7 ) f| Ut=o P 7~* ( J i)) and ^ stands for I- A. Let m fc = fi [\J k ~^ P~ { ( J 7 ) 
for k > 2. Then by mixing property ®, for sufficiently large ko G IN, 

Mfe « <5m fc <5m fco (1 - <5) fc ~ fc ° , fc > k . (9) 

The combination of (jHJ) and @ yields 

00 

iV 7 = S fco + m ko 5 ^(k + k ) (1 - 5) k ~ k ° , (10) 
fc=i 

where Sfc stands for the first ko terms on the right-hand side of (|5|>. By taking into account, 
S fco = O (1) and (fe + ko) (1 - <5) fc = O (S' 2 ) , from dJUJ) we obtain iV 7 = O ((T 1 ) . In a 

small neighborhood around the point of maximum, the graph of P 7 is to leading order quadratic. 

Therefore, the size of the window J 7 , 8 = fj, (J 7 ) = O (\/7~~7c) > an d ^7 = O ^(7 — 7 C )~ 2_ ^- To 
estimate iV~ in a larger neighborhood of 7 C , we need to review certain facts about the structure of P 7 . 
It follows from the construction of P 7 in T , that outside of the exponentially small neighborhood 

of WHC, 7, P 7 can be approximated by a linear map. This implies that for 7 > 7 = 7 C + O ( e_ " i ) 

the size of J 7 grows approximately linearly with 7, 5 ~ O (7 — 7) and iV" 7 = O ^(7 — 7) ^ for 
7 > 7. Therefore, in the vicinity of 7 C , there are two regions of qualitatively distinct asymptotic 
behaviors of iV 7 as a function of the distance from criticality. The numerical results shown in Figure 
0] confirm this conclusion and clearly indicate the boundary between these two regions, 7 ~ 2.943. 
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Figure 5: The histograms for the number of spikes within one burst and the mean values of the 
interspike intervals within one burst, T 7 . 

The proposed mechanism for irregular bursting implies that near the transition to bursting, iV 7 has 
a geometric distribution, whose parameters are determined by the width of the window of escape, 
J 7 (see Figure EJ). 

Next, we turn to estimating the expected value for the interspike intervals within one burst, T 7 . 
For small 7 — 7 C > 0, we have 

T 7 « / T(w, 7 c ) dn(w), (11) 

Jl-J-y 

where T (w, 7) stands for the period of oscillations in the fast subsystem Q,©- By the Lebesgue- 
Besicovitch differentiation theorem Jj T (w, 7 C ) dfj,(w) = O (fi ( J 7 )) = O (5), for small 7 — j c > 
0. Since 5 = O (-^7 — 7 C ) for 7 G (7 C , 7), in this region, we have T 7 = T 7 c — O (1/7 — 7 C ) . To 
estimate T 7 for 7 > 7, we again use the proximity of P 7 to a linear map in I — I, which implies 
that / — J 7 rs (u) , — O (^)) and dp,(w) ~ C2dw in I — I for some C2 > 0. In addition, by the 
well known results of the bifurcation theory ^ , 

T (to, 7 C ) T - C 3 log (wfjc - w) , 

where positive constants To and C3 are independent from w, because the fast subsystem is close to 
a homoclinic bifurcation. Using these approximations and (|11|). we obtain 

T y ^C 2 T(w,i c )dw^T -C 4 5\ln5\ + 0(5), 

J WO 

where positive constants T and C§ are independent of 5 and 5 = O (7 — 7 C ) . The numerical results in 
Figure \5\ show that the asymptotic behaviors of T 7 are different for 7 < 7 and 7 > 7. The sub linear 
character of T 7 in the latter region is consistent with our estimate T 7 — T = O ((7 — 7) |log (7 — 7)|) . 

The map-based approach to bursting employed in this paper reduces the problem of transition 
from (tonic) spiking to bursting to the analysis of bifurcation scenarios in the families of the first 
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Figure 6: The mechanism of the transition to bursting through a SN bifurcation. The SN bi- 
furcation gives rise to a periodic orbit corresponding to a bursting solution. The global return 
mechanism is schematically indicated by the dashed arc. 

return maps. The first event in these scenarios is the loss of stability of the fixed point corre- 
sponding to tonic spiking in the continuous system. For ID maps, there are two co-dimension 1 
bifurcations of fixed points: a SN and a PD bifurcation 4j. The bifurcation scenario described 
in [7] and in the present paper originates from a PD bifurcation. A complimentary mechanism of 
the transition to bursting is realized through a SN bifurcation (see Figure EJ). It has been recently 
demonstrated for a model for the leech heart interneuron ^3] (see also Remark 6.1(b) in 7 ). In 
analogy to the classification of excitability in spiking neuron models due to Rinzel and Ermentrout 
|10| . we propose to distinguish mechanisms of transition to bursting according to the underlying 
bifurcation mechanisms in the families of the first return maps: Type I transition is realized via 
a SN bifurcation of the fixed point (which coincides with the global bifurcation of a periodic orbit, 
see Figure El); Type II transition is realized through the disappearance (deflation) of the chaotic 
attractor (an attracting invariant interval with mixing property) preceded by a PD bifurcation. 
In both scenarios, the transition to bursting takes place after global bifurcations. However, the 
essential part of the bifurcation mechanism, which determines the traits of the bifurcating burst- 
ing patterns, in each case has a local character: a SN bifurcation in the Type I scenario and the 
deflation of the chaotic attractor in Type II. The bursting patterns arising through these mecha- 
nisms possess well-defined statistical properties in terms of the mean burst duration and interspike 
interval distributions within one burst. In particular, the interspike interval distributions in Type 
II bursting patterns are characterized by high variability, whereas those in Type I are localized. 
The mechanisms for transition to bursting in neuronal models are reminiscent to those studied in 
the context of the transition to turbulence via intermittency 0. In fact, the duration of (regu- 
lar) bursting near the transition in the Type I scenario and that of the laminar phases in the 
Type-I intermittency ^2] share the same asymptotics due to the proximity to a SN bifurcation in 
both cases. In contrast, the statistical properties of the bursting patterns in Type II scenario are 
determined by the long intervals of irregular (chaotic) behavior and have not been studied before. 

The author thanks to anonymous referees for insisting on careful numerical verification of the 
analytical estimates presented in this paper. This work was partially supported by the National 
Science Foundation under Grant No. 0417624. 
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